rm(list = ls())
library(pacman)
p_load(lme4,lmerTest ,performance, tidyverse,lubridate, lattice, devtools, VGAM, ggpubr, plotly)
library("knitr")
name <- data0.anal %>% distinct(iso_code)
# JOR is inconsistent data
hist(data0.anal$week.case.eff)
check <- data0.anal %>% filter(iso_code == "JOR") %>% select(iso_code, time, iso_code, total.cases, week.case.eff)
summary(data0.anal$school.close.red)
## 0 2 3
## 68 61 471
summary(data0.anal$work.close.red)
## 0 2 3
## 156 313 131
summary(data0.anal$restric.gather.red)
## 0 1 3 4
## 139 73 53 335
setwd("F:/STUDY HCA IN TAIWAN/0. ideas for thesis/COVID19 and NPIs/Data analysis/dataset/data_sep")
week20 <- data0.anal %>% filter(time == "19") %>% select(country, population, total.cases)
## Adding missing grouping variables: `iso_code`
write.table(week20, file = "tt_case_2020_w20.csv", sep = ",")
check.time <- data0.anal %>% select(id, week.case.eff, growth.rate.manual, timeStart, timeEff)
## Adding missing grouping variables: `iso_code`
# Group plot
plot.group.case <- ggplot(data0.anal, aes(time+2, week.case.eff)) +
stat_smooth(method = "loess", se = F, size = 2) +
stat_summary(aes(col = country), fun.y = mean, geom = "line", alpha = 0.5) +
xlab("Week") +
ylab(" Average weekly growth rate")+
theme_bw() +
theme(strip.background = element_blank()) +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
theme(legend.position = "bottom") +
scale_x_continuous(breaks = seq(2, 25, by = 1))
## Warning: `fun.y` is deprecated. Use `fun` instead.
plot.group.case
## `geom_smooth()` using formula 'y ~ x'
plot.group.string <- ggplot(data0.anal, aes(time, string.idx)) +
stat_smooth(method = "loess", se = F, size = 2) +
stat_summary(aes(col = country), fun.y = mean, geom = "line", alpha = 0.5) +
xlab("Week") +
ylab(" Stringency index")+
theme_bw() +
theme(strip.background = element_blank()) +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
theme(legend.position = "bottom") +
scale_x_continuous(breaks = seq(0, 23, by = 1))
## Warning: `fun.y` is deprecated. Use `fun` instead.
plot.group.string
## `geom_smooth()` using formula 'y ~ x'
cairo_pdf(filename = "F:/STUDY HCA IN TAIWAN/0. ideas for thesis/COVID19 and NPIs/Data analysis/results/plot_sep/fig1_NPI2020_Sep29.pdf",
width = 18, height = 8, #inch
onefile = TRUE, family = "Segoe UI")
merge1 <- ggarrange(plot.group.string, plot.group.case, hjust = -1,
ncol = 2, nrow = 1)
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
merge1
dev.off()
## png
## 2
ggplotly(plot.group.case)
## `geom_smooth()` using formula 'y ~ x'
ggplotly(plot.group.string)
## `geom_smooth()` using formula 'y ~ x'
#https://www.guru99.com/r-scatter-plot-ggplot2.html#2
# clean daily grow rate = 0, and NAs
## NAs = 36 obs
## zero growth rate = 58
## Check zero growth rate
data0.anal1 <- data0.anal
check <- data0.anal1 %>% filter(week.case.eff == 0 )
check1 <- data0.anal1 %>% filter(week.case.eff != 0 ) %>% select(id, week.case.eff)
## Adding missing grouping variables: `iso_code`
## Check NAs
summary(data0.anal$week.case.eff)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000000 0.004217 0.017600 0.039388 0.051100 0.609670
#probit transformation
data0.anal <- data0.anal %>% filter(week.case.eff != 0) %>% mutate(probit.case = probitlink( week.case.eff))
qqnorm(data0.anal$probit.case)
hist(data0.anal$probit.case)
boxplot(data0.anal$probit.case)
summary(data0.anal$test.policy.red)
## Warning: Unknown or uninitialised column: `test.policy.red`.
## Length Class Mode
## 0 NULL NULL
data0.anal <- data0.anal %>% mutate(test.policy.red = case_when(
test.policy==0~0,
test.policy==1~0,
test.policy==2~1,
test.policy==3~2
))
data0.anal$test.policy.red <- factor(data0.anal$test.policy.red)
summary(data0.anal$inter.travel)
## 0 1 2 3 4
## 22 9 37 202 270
data0.anal <- data0.anal %>% mutate(inter.travel.red = case_when(
inter.travel==0~0,
inter.travel==1~0,
inter.travel==2~0,
inter.travel==3~1,
inter.travel==4~2
))
data0.anal$inter.travel.red <- factor(data0.anal$inter.travel.red)
summary(data0.anal$wear.mask)
## 0 1 2 3 4
## 188 58 101 90 103
data0.anal <- data0.anal %>% mutate(wear.mask.red = case_when(
wear.mask==0~0,
wear.mask==1~1,
wear.mask==2~1,
wear.mask==3~2,
wear.mask==4~3
))
data0.anal$wear.mask.red <- factor(data0.anal$wear.mask.red)
data0.anal <- data0.anal %>% mutate(wear.mask.red1 = case_when(
wear.mask==0~0,
wear.mask==1~0,
wear.mask==2~1,
wear.mask==3~2,
wear.mask==4~3
))
data0.anal$wear.mask.red1 <- factor(data0.anal$wear.mask.red1)
data0.anal$pop_log <- log(data0.anal$population)
hist(data0.anal$pop_log)
data0.anal$pop_denlog <- log(data0.anal$pop_density)
hist(data0.anal$pop_denlog)
data0.anal$median_age.log <- log(data0.anal$median_age)
hist(data0.anal$median_age.log)
data0.anal$gdp.log <- log(data0.anal$gdp)
hist(data0.anal$gdp.log)
data0.anal$incomegroup1 <- factor(data0.anal$incomegroup, levels = c("Low income", "Lower middle income", "Upper middle income" , "High income"))
data0.anal <- data0.anal %>% mutate(incomegroup2 = case_when(
incomegroup1 == "Low income" ~ "LMIC",
incomegroup1 == "Lower middle income" ~ "LMIC",
incomegroup1 == "Upper middle income" ~ "upper_income",
incomegroup1 == "High income" ~ "high_income"
))
data0.anal$incomegroup2 <- factor(data0.anal$incomegroup2, levels = c( "LMIC", "upper_income","high_income"))
library(performance)
aic <- function(model){
fix2 <- performance(model)
x<- as.character(substitute(model))
fix2 <- fix2 %>%
mutate( npi = x)
}
model.school.cls <- lmer(probit.case ~ school.close.red +
time + (time|iso_code),
data= data0.anal, REML = F, lmerControl(optimizer = "bobyqa") )
summary(model.school.cls)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: probit.case ~ school.close.red + time + (time | iso_code)
## Data: data0.anal
## Control: lmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 484.9 519.3 -234.5 468.9 532
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.1656 -0.5395 -0.0463 0.4680 4.5459
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## iso_code (Intercept) 0.2044411 0.45215
## time 0.0006648 0.02578 -0.42
## Residual 0.1066855 0.32663
## Number of obs: 540, groups: iso_code, 30
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -1.230536 0.106689 57.448375 -11.534 < 2e-16 ***
## school.close.red2 0.054009 0.094487 517.578250 0.572 0.56784
## school.close.red3 -0.211454 0.069518 497.657722 -3.042 0.00248 **
## time -0.063180 0.005503 30.598704 -11.480 1.28e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) sch..2 sch..3
## schl.cls.r2 -0.344
## schl.cls.r3 -0.558 0.662
## time -0.385 -0.171 -0.035
x1 <- aic(model.school.cls)
model.work.cls <- lmer(probit.case ~ work.close.red +
time + (time|iso_code),
data= data0.anal, REML = F, lmerControl(optimizer = "bobyqa") )
summary(model.work.cls)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: probit.case ~ work.close.red + time + (time | iso_code)
## Data: data0.anal
## Control: lmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 488.9 523.3 -236.5 472.9 532
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.1447 -0.5262 -0.0373 0.4362 4.5904
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## iso_code (Intercept) 0.2276444 0.47712
## time 0.0007713 0.02777 -0.53
## Residual 0.1070298 0.32715
## Number of obs: 540, groups: iso_code, 30
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -1.299756 0.098404 35.768193 -13.208 2.53e-15 ***
## work.close.red2 -0.156166 0.057166 399.519932 -2.732 0.00658 **
## work.close.red3 -0.215738 0.053792 509.475989 -4.011 6.96e-05 ***
## time -0.059369 0.005857 30.094125 -10.136 3.23e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) wrk..2 wrk..3
## wrk.cls.rd2 -0.287
## wrk.cls.rd3 -0.310 0.656
## time -0.502 -0.173 -0.023
x2<- aic(model.work.cls)
model.pubevent <- lmer(probit.case ~ cancel.pubevent.red +
time + (time|iso_code),
data= data0.anal, REML = F, lmerControl(optimizer = "bobyqa") )
summary(model.pubevent)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: probit.case ~ cancel.pubevent.red + time + (time | iso_code)
## Data: data0.anal
## Control: lmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 494.6 524.7 -240.3 480.6 533
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.9901 -0.5414 -0.0326 0.4599 4.3079
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## iso_code (Intercept) 0.2179870 0.46689
## time 0.0006557 0.02561 -0.48
## Residual 0.1094033 0.33076
## Number of obs: 540, groups: iso_code, 30
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -1.302306 0.100308 40.951960 -12.983 4.19e-16 ***
## cancel.pubevent.red2 -0.163539 0.056440 477.439874 -2.898 0.00393 **
## time -0.059070 0.005393 28.857875 -10.953 8.56e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) cnc..2
## cncl.pbvn.2 -0.417
## time -0.461 -0.064
x3<- aic(model.pubevent)
model.restr.gather <- lmer(probit.case ~ restric.gather.red +
time + (time|iso_code),
data= data0.anal, REML = F, lmerControl(optimizer = "bobyqa") )
summary(model.restr.gather)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: probit.case ~ restric.gather.red + time + (time | iso_code)
## Data: data0.anal
## Control: lmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 492.6 531.3 -237.3 474.6 531
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.9120 -0.5585 -0.0522 0.4664 4.6081
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## iso_code (Intercept) 0.2244701 0.47378
## time 0.0006785 0.02605 -0.44
## Residual 0.1074566 0.32781
## Number of obs: 540, groups: iso_code, 30
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -1.302949 0.098989 37.681114 -13.163 1.14e-15 ***
## restric.gather.red1 -0.035307 0.083952 524.739493 -0.421 0.674243
## restric.gather.red3 -0.147442 0.088206 500.481696 -1.672 0.095234 .
## restric.gather.red4 -0.222364 0.060163 522.223255 -3.696 0.000242 ***
## time -0.057899 0.005508 29.903788 -10.512 1.46e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) rst..1 rst..3 rst..4
## rstrc.gth.1 -0.204
## rstrc.gth.3 -0.238 0.406
## rstrc.gth.4 -0.353 0.384 0.554
## time -0.417 -0.108 -0.128 -0.128
x4 <- aic(model.restr.gather)
model.cls.trans <- lmer(probit.case ~ close.transport.red +
time + (time|iso_code),
data= data0.anal, REML = F, lmerControl(optimizer = "bobyqa") )
summary(model.cls.trans)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: probit.case ~ close.transport.red + time + (time | iso_code)
## Data: data0.anal
## Control: lmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 499.0 529.0 -242.5 485.0 533
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.8582 -0.5388 -0.0419 0.4456 4.6640
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## iso_code (Intercept) 0.2474150 0.49741
## time 0.0007383 0.02717 -0.55
## Residual 0.1096458 0.33113
## Number of obs: 540, groups: iso_code, 30
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -1.368495 0.100817 34.658429 -13.57 1.96e-15 ***
## close.transport.red1 -0.093965 0.049971 514.440037 -1.88 0.0606 .
## time -0.059710 0.005648 28.885222 -10.57 1.93e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) cls..1
## cls.trnsp.1 -0.288
## time -0.559 -0.038
x5 <- aic(model.cls.trans)
model.stay.home <- lmer(probit.case ~ stay.home.red +
time + (time|iso_code),
data= data0.anal, REML = F, lmerControl(optimizer = "bobyqa") )
summary(model.stay.home)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: probit.case ~ stay.home.red + time + (time | iso_code)
## Data: data0.anal
## Control: lmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 492.5 522.5 -239.2 478.5 533
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.0139 -0.5302 -0.0375 0.4581 4.5679
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## iso_code (Intercept) 0.2459927 0.49598
## time 0.0007357 0.02712 -0.51
## Residual 0.1078785 0.32845
## Number of obs: 540, groups: iso_code, 30
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -1.366991 0.097769 31.340528 -13.982 5.13e-15 ***
## stay.home.red1 -0.150050 0.046519 533.691032 -3.226 0.00133 **
## time -0.059139 0.005633 29.064622 -10.498 2.11e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) sty..1
## stay.hm.rd1 -0.178
## time -0.543 -0.053
x6 <- aic(model.stay.home)
model.restr.move <- lmer(probit.case ~ restric.move.red +
time + (time|iso_code),
data= data0.anal, REML = F, lmerControl(optimizer = "bobyqa") )
summary(model.restr.move)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: probit.case ~ restric.move.red + time + (time | iso_code)
## Data: data0.anal
## Control: lmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 497.0 527.0 -241.5 483.0 533
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.0133 -0.5325 -0.0323 0.4478 4.5814
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## iso_code (Intercept) 0.2481191 0.49812
## time 0.0007701 0.02775 -0.57
## Residual 0.1092041 0.33046
## Number of obs: 540, groups: iso_code, 30
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -1.371390 0.099162 32.374319 -13.83 3.88e-15 ***
## restric.move.red2 -0.103169 0.043910 532.547488 -2.35 0.0192 *
## time -0.059713 0.005741 28.756405 -10.40 2.95e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) rst..2
## rstrc.mv.r2 -0.222
## time -0.592 -0.030
x7 <- aic(model.restr.move)
model.inter.trav <- lmer(probit.case ~ inter.travel.red +
time + (time|iso_code),
data= data0.anal, REML = F, lmerControl(optimizer = "bobyqa") )
summary(model.inter.trav)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: probit.case ~ inter.travel.red + time + (time | iso_code)
## Data: data0.anal
## Control: lmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 482.4 516.7 -233.2 466.4 532
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.8349 -0.5442 -0.0098 0.4428 4.7752
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## iso_code (Intercept) 0.2593476 0.50926
## time 0.0006816 0.02611 -0.46
## Residual 0.1048992 0.32388
## Number of obs: 540, groups: iso_code, 30
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -1.284909 0.112827 49.711369 -11.388 1.84e-15 ***
## inter.travel.red1 -0.050367 0.071978 497.830959 -0.700 0.484413
## inter.travel.red2 -0.286159 0.075038 485.356616 -3.814 0.000155 ***
## time -0.058036 0.005462 29.082018 -10.626 1.58e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) int..1 int..2
## intr.trvl.1 -0.453
## intr.trvl.2 -0.458 0.727
## time -0.419 -0.040 -0.080
x8 <- aic(model.inter.trav)
model.test <- lmer(probit.case ~ test.policy.red +
time + (time|iso_code),
data= data0.anal, REML = F, lmerControl(optimizer = "bobyqa") )
summary(model.test)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: probit.case ~ test.policy.red + time + (time | iso_code)
## Data: data0.anal
## Control: lmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 495.2 529.5 -239.6 479.2 532
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.8244 -0.5376 -0.0278 0.4330 4.7334
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## iso_code (Intercept) 0.248623 0.49862
## time 0.000638 0.02526 -0.61
## Residual 0.109687 0.33119
## Number of obs: 540, groups: iso_code, 30
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -1.38791 0.09984 32.96732 -13.901 2.42e-15 ***
## test.policy.red1 -0.08567 0.06877 411.79676 -1.246 0.21355
## test.policy.red2 -0.25225 0.08171 385.83663 -3.087 0.00217 **
## time -0.05424 0.00572 36.55832 -9.481 2.17e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) tst..1 tst..2
## tst.plcy.r1 -0.247
## tst.plcy.r2 -0.132 0.486
## time -0.511 -0.261 -0.345
x9 <- aic(model.test)
model.tracing <- lmer(probit.case ~ contact.tracing +
time + (time|iso_code),
data= data0.anal, REML = F, lmerControl(optimizer = "bobyqa") )
summary(model.tracing)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: probit.case ~ contact.tracing + time + (time | iso_code)
## Data: data0.anal
## Control: lmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 493.2 527.6 -238.6 477.2 532
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.8719 -0.5219 -0.0381 0.4368 4.7427
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## iso_code (Intercept) 0.3193535 0.56511
## time 0.0007502 0.02739 -0.70
## Residual 0.1080664 0.32873
## Number of obs: 540, groups: iso_code, 30
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -1.806154 0.252055 128.554319 -7.166 5.37e-11 ***
## contact.tracing1 0.591454 0.251570 210.954308 2.351 0.0196 *
## contact.tracing2 0.293960 0.236264 187.365946 1.244 0.2150
## time -0.060189 0.005721 29.197895 -10.521 1.90e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) cntc.1 cntc.2
## cntct.trcn1 -0.880
## cntct.trcn2 -0.894 0.934
## time -0.211 -0.087 -0.109
x10 <- aic(model.tracing)
model.mask <- lmer(probit.case ~ wear.mask.red +
time + (time|iso_code),
data= data0.anal, REML = F, lmerControl(optimizer = "bobyqa") )
summary(model.mask)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: probit.case ~ wear.mask.red + time + (time | iso_code)
## Data: data0.anal
## Control: lmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 482.8 521.4 -232.4 464.8 531
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.3769 -0.5310 -0.0397 0.4861 4.4528
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## iso_code (Intercept) 0.2190665 0.46805
## time 0.0006316 0.02513 -0.65
## Residual 0.1078840 0.32846
## Number of obs: 540, groups: iso_code, 30
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -1.356786 0.092316 29.572995 -14.697 3.94e-15 ***
## wear.mask.red1 -0.343768 0.068662 493.118928 -5.007 7.72e-07 ***
## wear.mask.red2 -0.223732 0.085326 332.601163 -2.622 0.00914 **
## wear.mask.red3 -0.125477 0.071999 404.955827 -1.743 0.08214 .
## time -0.050290 0.005932 42.460894 -8.478 1.11e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) wr.m.1 wr.m.2 wr.m.3
## wer.msk.rd1 -0.143
## wer.msk.rd2 -0.065 0.446
## wer.msk.rd3 -0.042 0.285 0.300
## time -0.558 -0.299 -0.354 -0.344
x11 <- aic(model.mask)
model.info <- lmer(probit.case ~ info.campaign +
time + (time|iso_code),
data= data0.anal, REML = F, lmerControl(optimizer = "bobyqa") )
summary(model.info)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: probit.case ~ info.campaign + time + (time | iso_code)
## Data: data0.anal
## Control: lmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 500.9 535.2 -242.5 484.9 532
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.8260 -0.5417 -0.0377 0.4489 4.6673
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## iso_code (Intercept) 0.2448265 0.49480
## time 0.0006834 0.02614 -0.56
## Residual 0.1101943 0.33196
## Number of obs: 540, groups: iso_code, 30
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -1.609405 0.217299 280.481998 -7.406 1.53e-12 ***
## info.campaign1 0.428322 0.251262 495.424675 1.705 0.0889 .
## info.campaign2 0.172464 0.201308 462.073047 0.857 0.3920
## time -0.059719 0.005522 29.434366 -10.815 9.03e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) inf.c1 inf.c2
## info.cmpgn1 -0.745
## info.cmpgn2 -0.895 0.795
## time -0.191 -0.017 -0.088
x12 <- aic(model.info)
#population: no
model.pop <- lmer(probit.case ~ pop_log +
time + (time|iso_code),
data= data0.anal, REML = F, lmerControl(optimizer = "bobyqa") )
summary(model.pop)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: probit.case ~ pop_log + time + (time | iso_code)
## Data: data0.anal
## Control: lmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 501.3 531.4 -243.7 487.3 533
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.8339 -0.5429 -0.0345 0.4384 4.6547
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## iso_code (Intercept) 0.2516624 0.50166
## time 0.0007352 0.02711 -0.60
## Residual 0.1106882 0.33270
## Number of obs: 540, groups: iso_code, 30
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -0.686257 0.698982 31.227979 -0.982 0.334
## pop_log -0.042756 0.040183 30.539823 -1.064 0.296
## time -0.060179 0.005643 28.806722 -10.664 1.63e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) pop_lg
## pop_log -0.990
## time -0.097 0.008
y1 <- aic(model.pop)
#population density
model.popden <- lmer(probit.case ~ pop_denlog +
time + (time|iso_code),
data= data0.anal, REML = F, lmerControl(optimizer = "bobyqa") )
summary(model.popden)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: probit.case ~ pop_denlog + time + (time | iso_code)
## Data: data0.anal
## Control: lmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 501.5 531.6 -243.8 487.5 533
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.8336 -0.5381 -0.0393 0.4384 4.6406
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## iso_code (Intercept) 0.2391485 0.48903
## time 0.0007352 0.02711 -0.57
## Residual 0.1106907 0.33270
## Number of obs: 540, groups: iso_code, 30
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -1.652227 0.252408 33.689187 -6.546 1.77e-07 ***
## pop_denlog 0.045603 0.046536 30.189707 0.980 0.335
## time -0.060053 0.005641 28.813550 -10.646 1.69e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) pp_dnl
## pop_denlog -0.926
## time -0.240 0.009
y2 <- aic(model.popden)
# median age
model.age <- lmer(probit.case ~ median_age +
time + (time|iso_code),
data= data0.anal, REML = F, lmerControl(optimizer = "bobyqa") )
summary(model.age)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: probit.case ~ median_age + time + (time | iso_code)
## Data: data0.anal
## Control: lmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 495.7 525.8 -240.9 481.7 533
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.8248 -0.5395 -0.0350 0.4409 4.6803
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## iso_code (Intercept) 0.2602538 0.5102
## time 0.0007343 0.0271 -0.71
## Residual 0.1107129 0.3327
## Number of obs: 540, groups: iso_code, 30
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -0.477093 0.341241 34.144716 -1.398 0.17110
## median_age -0.030262 0.010459 30.882357 -2.893 0.00693 **
## time -0.060290 0.005645 28.726092 -10.681 1.62e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) medn_g
## median_age -0.957
## time -0.217 0.007
y3 <- aic(model.age)
# gdp
model.gdp <- lmer(probit.case ~ gdp.log +
time + (time|iso_code),
data= data0.anal, REML = F, lmerControl(optimizer = "bobyqa") )
summary(model.gdp)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: probit.case ~ gdp.log + time + (time | iso_code)
## Data: data0.anal
## Control: lmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 502 532 -244 488 533
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.8315 -0.5445 -0.0348 0.4372 4.6586
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## iso_code (Intercept) 0.2587356 0.5087
## time 0.0007347 0.0271 -0.61
## Residual 0.1106916 0.3327
## Number of obs: 540, groups: iso_code, 30
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -0.910145 0.733937 31.056382 -1.240 0.224
## gdp.log -0.053077 0.075327 30.370869 -0.705 0.486
## time -0.060187 0.005642 28.790591 -10.667 1.63e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) gdp.lg
## gdp.log -0.991
## time -0.098 0.012
y4 <- aic(model.gdp)
# test per 1000 => too many missing because the early stage of pandemic
model.test1000 <- lmer(probit.case ~ test.thousand +
time + (time|iso_code),
data= data0.anal, REML = F, lmerControl(optimizer = "bobyqa") )
summary(model.test1000)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: probit.case ~ test.thousand + time + (time | iso_code)
## Data: data0.anal
## Control: lmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 301.3 329.6 -143.7 287.3 411
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.3822 -0.4110 -0.0188 0.3532 5.1613
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## iso_code (Intercept) 0.1956395 0.44231
## time 0.0009157 0.03026 -0.64
## Residual 0.0879862 0.29662
## Number of obs: 418, groups: iso_code, 24
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -1.3520967 0.0986303 23.2956823 -13.709 1.21e-12 ***
## test.thousand -0.0009576 0.0008348 44.7984633 -1.147 0.257
## time -0.0577606 0.0075758 26.9090064 -7.624 3.43e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) tst.th
## test.thosnd 0.099
## time -0.665 -0.414
y5 <- aic(model.test1000)
# percentage of urban
model.per_urban <- lmer(probit.case ~ per_urban +
time + (time|iso_code),
data= data0.anal, REML = F, lmerControl(optimizer = "bobyqa") )
summary(model.per_urban)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: probit.case ~ per_urban + time + (time | iso_code)
## Data: data0.anal
## Control: lmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 502.0 532.1 -244.0 488.0 533
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.8322 -0.5443 -0.0350 0.4377 4.6549
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## iso_code (Intercept) 0.2530655 0.50306
## time 0.0007347 0.02711 -0.59
## Residual 0.1106905 0.33270
## Number of obs: 540, groups: iso_code, 30
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -1.307725 0.204141 36.795138 -6.406 1.82e-07 ***
## per_urban -0.001889 0.002944 30.543970 -0.642 0.526
## time -0.060181 0.005642 28.816358 -10.667 1.61e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) pr_rbn
## per_urban -0.878
## time -0.315 0.015
y6 <- aic(model.per_urban)
# uhc_index
model.uhc <- lmer(probit.case ~ uhc_index +
time + (time|iso_code),
data= data0.anal, REML = F, lmerControl(optimizer = "bobyqa") )
summary(model.uhc)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: probit.case ~ uhc_index + time + (time | iso_code)
## Data: data0.anal
## Control: lmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 494.8 524.9 -240.4 480.8 533
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.8284 -0.5384 -0.0308 0.4440 4.6733
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## iso_code (Intercept) 0.2347485 0.48451
## time 0.0007333 0.02708 -0.68
## Residual 0.1107152 0.33274
## Number of obs: 540, groups: iso_code, 30
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -0.244745 0.404956 33.184394 -0.604 0.54970
## uhc_index -0.017419 0.005823 30.926225 -2.992 0.00541 **
## time -0.060373 0.005639 28.815358 -10.706 1.48e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) uhc_nd
## uhc_index -0.972
## time -0.178 0.015
y7 <- aic(model.uhc)
# health worker density
model.hworker <- lmer(probit.case ~ healthworker_den +
time + (time|iso_code),
data= data0.anal, REML = F, lmerControl(optimizer = "bobyqa") )
summary(model.hworker)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: probit.case ~ healthworker_den + time + (time | iso_code)
## Data: data0.anal
## Control: lmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 499.2 529.3 -242.6 485.2 533
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.8323 -0.5420 -0.0319 0.4387 4.6618
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## iso_code (Intercept) 0.2443999 0.49437
## time 0.0007337 0.02709 -0.63
## Residual 0.1107056 0.33272
## Number of obs: 540, groups: iso_code, 30
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -1.166817 0.168425 40.689289 -6.928 2.17e-08 ***
## healthworker_den -0.004002 0.002162 31.275638 -1.851 0.0737 .
## time -0.060325 0.005640 28.835170 -10.697 1.50e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) hlthw_
## hlthwrkr_dn -0.821
## time -0.392 0.019
y8 <- aic(model.hworker)
#Income group
model.income <- lmer(probit.case ~ incomegroup2 +
time + (time|iso_code),
data= data0.anal, REML = F, lmerControl(optimizer = "bobyqa") )
summary(model.income)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: probit.case ~ incomegroup2 + time + (time | iso_code)
## Data: data0.anal
## Control: lmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 500.2 534.6 -242.1 484.2 532
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.8123 -0.5392 -0.0420 0.4475 4.6250
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## iso_code (Intercept) 0.2109792 0.45932
## time 0.0007316 0.02705 -0.57
## Residual 0.1107088 0.33273
## Number of obs: 540, groups: iso_code, 30
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -1.297804 0.123201 38.978726 -10.534 5.78e-13 ***
## incomegroup2upper_income -0.352801 0.168818 30.379305 -2.090 0.0451 *
## incomegroup2high_income -0.061070 0.174839 30.322594 -0.349 0.7293
## time -0.060232 0.005628 28.965326 -10.703 1.40e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) incmgrp2p_ incmgrp2h_
## incmgrp2pp_ -0.582
## incmgrp2hg_ -0.562 0.405
## time -0.462 0.018 0.017
y9 <- aic(model.income)
# SDI
model.sdi <- lmer(probit.case ~ sdi_index +
time + (time|iso_code),
data= data0.anal, REML = F, lmerControl(optimizer = "bobyqa") )
summary(model.sdi)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: probit.case ~ sdi_index + time + (time | iso_code)
## Data: data0.anal
## Control: lmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 501.5 531.6 -243.8 487.5 533
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.8302 -0.5449 -0.0345 0.4348 4.6649
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## iso_code (Intercept) 0.2622524 0.5121
## time 0.0007344 0.0271 -0.62
## Residual 0.1106941 0.3327
## Number of obs: 540, groups: iso_code, 30
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -1.060386 0.371535 33.263400 -2.854 0.00737 **
## sdi_index -0.542060 0.536009 30.775256 -1.011 0.31976
## time -0.060230 0.005642 28.779814 -10.674 1.61e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) sd_ndx
## sdi_index -0.964
## time -0.188 0.014
y10 <- aic(model.sdi)
# mobi.composite
model.mobi <- lmer(probit.case ~ mobi.composite +
time + (time|iso_code),
data= data0.anal, REML = F, lmerControl(optimizer = "bobyqa") )
summary(model.mobi)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: probit.case ~ mobi.composite + time + (time | iso_code)
## Data: data0.anal
## Control: lmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 393.6 422.9 -189.8 379.6 478
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.2626 -0.5368 -0.0019 0.4520 4.7071
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## iso_code (Intercept) 0.2106074 0.45892
## time 0.0005974 0.02444 -0.29
## Residual 0.0977931 0.31272
## Number of obs: 485, groups: iso_code, 25
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -1.235e+00 1.024e-01 3.117e+01 -12.064 2.81e-13 ***
## mobi.composite 3.613e-03 8.774e-04 4.519e+02 4.118 4.54e-05 ***
## time -6.415e-02 5.528e-03 2.454e+01 -11.604 1.88e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) mb.cmp
## mobi.compst 0.341
## time -0.369 -0.072
y11 <- aic(model.mobi)
aic_table <- rbind(x1, x2, x3, x4, x6, x7, x8, x9, x10, x11, y3, y7)
model.mul <- lmer(probit.case ~ uhc_index + inter.travel.red + wear.mask.red + test.policy.red +
time + (time|iso_code),
data= data0.anal, REML = F, lmerControl(optimizer = "bobyqa") )
summary(model.mul)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: probit.case ~ uhc_index + inter.travel.red + wear.mask.red +
## test.policy.red + time + (time | iso_code)
## Data: data0.anal
## Control: lmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 464.5 524.6 -218.3 436.5 526
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.3417 -0.5247 -0.0151 0.4634 4.6497
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## iso_code (Intercept) 0.2186297 0.46758
## time 0.0005131 0.02265 -0.64
## Residual 0.1025931 0.32030
## Number of obs: 540, groups: iso_code, 30
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -0.319755 0.423282 33.223032 -0.755 0.45532
## uhc_index -0.013109 0.006120 31.463177 -2.142 0.04003 *
## inter.travel.red1 -0.059617 0.071211 494.020215 -0.837 0.40289
## inter.travel.red2 -0.238736 0.074238 461.469385 -3.216 0.00139 **
## wear.mask.red1 -0.311671 0.068913 506.738186 -4.523 7.61e-06 ***
## wear.mask.red2 -0.184827 0.084831 333.205558 -2.179 0.03005 *
## wear.mask.red3 -0.089169 0.070741 335.209161 -1.260 0.20837
## test.policy.red1 -0.092283 0.067384 374.048339 -1.370 0.17166
## test.policy.red2 -0.169078 0.081638 337.485494 -2.071 0.03911 *
## time -0.046106 0.005828 52.223350 -7.912 1.71e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) uhc_nd int..1 int..2 wr.m.1 wr.m.2 wr.m.3 tst..1 tst..2
## uhc_index -0.966
## intr.trvl.1 -0.068 -0.055
## intr.trvl.2 -0.128 0.013 0.732
## wer.msk.rd1 0.073 -0.121 0.084 -0.003
## wer.msk.rd2 0.076 -0.113 0.159 0.126 0.470
## wer.msk.rd3 -0.068 0.059 0.053 0.005 0.288 0.313
## tst.plcy.r1 0.008 -0.064 0.001 -0.100 0.163 0.017 -0.040
## tst.plcy.r2 0.162 -0.175 -0.125 -0.196 0.000 -0.122 -0.151 0.502
## time -0.175 0.091 -0.063 -0.036 -0.327 -0.341 -0.299 -0.266 -0.263
performance(model.mul)
## # Indices of model performance
##
## AIC | BIC | R2 (cond.) | R2 (marg.) | ICC | RMSE | Sigma
## -------------------------------------------------------------------
## 464.528 | 524.610 | 0.761 | 0.405 | 0.598 | 0.305 | 0.320
# model.mul1 <- lmer(probit.case ~
# work.close.red + uhc_index + wear.mask + test.policy + incomegroup2 +
# time + (time|iso_code),
# data= data0.anal, REML = F, lmerControl(optimizer = "bobyqa") )
# summary(model.mul1)
# performance(model.mul1)
# anova(model.mul, model.mul1)
p_load(broom.mixed)
sum_table <- tidy(model.mul, conf.int = T)
# create new table summary of model
new_table <- sum_table[2:15,]
term1 <- c(
"UHC Index\n",
"Stay at home\nwith minimal exception\n",
"Screening arrivals\n",
"Quarantine arrivals\nfrom some or all regions\n",
"Ban arrivals from\nsome regions\n",
"Ban on all regions\n",
"School close some \nlevels\n",
"School close all levels\n",
"Only those who both have\nsymptoms+ meet specific criteria\n",
"Testing of anyone \nshowing symptoms\n",
"Open public testing\n",
"Upper income countries\n",
"High income countries\n",
"Time\n"
)
new_table1 <- cbind(new_table, term1)
new_table1$term2 <- factor(new_table1$term1, levels = new_table1$term1)
cairo_pdf(filename = "F:/STUDY HCA IN TAIWAN/0. ideas for thesis/COVID19 and NPIs/Data analysis/results/plot_sep/fig2_year2020_eff_NPI_Sep25.pdf",
width = 15, height = 10, #inch
onefile = TRUE, family = "Segoe UI")
p.year2020<- ggplot(new_table1, aes(x = term2, y = estimate,
ymin = conf.low, ymax = conf.high)) +
geom_hline( yintercept = 0, color = 'Blue' ) + ylim(-1.5,1) +
geom_linerange() + geom_point() + coord_flip() +
ylab("Probit Transformation of Average Weekly Growth Rate ") +
xlab("Factors") +
theme_minimal() +
theme(strip.background = element_rect(NA),
axis.text = element_text(size = 20),
axis.title = element_text(size = 20),
legend.text = element_text(size = 10),
legend.title = element_text(size = 10),
strip.text = element_text(size = 8),
legend.position = "top")
p.year2020
## Warning: Removed 5 rows containing missing values (geom_segment).
## Warning: Removed 1 rows containing missing values (geom_point).
dev.off()
## png
## 2
#devtools::install_github("goodekat/redres")
library(redres)
data0.anal1 <- data0.anal %>% filter( is.na(uhc_index) ==F)
rc_resids <- compute_redres(model.mul)
pm_resids <- compute_redres(model.mul, type = "pearson_mar")
sc_resids <- compute_redres(model.mul, type = "std_cond")
resids <- data.frame(data0.anal1$iso_code, rc_resids, pm_resids, sc_resids)
plot_redres(model.mul, type = "std_cond")
plot_resqq(model.mul)
plot_ranef(model.mul)
plot(data0.anal1$time, rc_resids,ylim=c(-3,3))
abline(h=0, col="blue")
plot(data0.anal1$iso_code, sc_resids)
abline(h=c(0,2.96,-2.96), col="blue")
# # run other diagnosis
# p_load(patchwork, randomForest)
#
# check_model(model.mul)
#
plot(check_distribution(model.mul))
## Warning: Removed 4 rows containing missing values (geom_segment).
## Warning: Removed 4 rows containing missing values (geom_point).
#
check_collinearity(model.mul)
## # Check for Multicollinearity
##
## Low Correlation
##
## Term VIF Increased SE Tolerance
## uhc_index 1.08 1.04 0.93
## inter.travel.red 1.10 1.05 0.91
## wear.mask.red 1.47 1.21 0.68
## test.policy.red 1.35 1.16 0.74
## time 1.43 1.20 0.70
#
pp_check(model.mul, 250)
library(glmmTMB)
## Warning in checkMatrixPackageVersion(): Package version inconsistency detected.
## TMB was built with Matrix version 1.3.4
## Current Matrix version is 1.3.3
## Please re-install 'TMB' from source using install.packages('TMB', type = 'source') or ask CRAN for a binary version of 'TMB' matching CRAN's 'Matrix' package
##
## Attaching package: 'glmmTMB'
## The following object is masked from 'package:VGAM':
##
## betabinomial
model.glm <- glmmTMB(week.case.eff ~ uhc_index + inter.travel.red + wear.mask.red + test.policy.red +
time + (time|iso_code), data = data0.anal, family = beta_family(link="probit"), REML = TRUE)
summary(model.glm)
## Family: beta ( probit )
## Formula:
## week.case.eff ~ uhc_index + inter.travel.red + wear.mask.red +
## test.policy.red + time + (time | iso_code)
## Data: data0.anal
##
## AIC BIC logLik deviance df.resid
## -2679.9 -2619.9 1354.0 -2707.9 536
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev. Corr
## iso_code (Intercept) 0.1224981 0.3500
## time 0.0002592 0.0161 -0.79
## Number of obs: 540, groups: iso_code, 30
##
## Dispersion parameter for beta family (): 43.3
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.654016 0.282376 -2.316 0.020552 *
## uhc_index -0.008340 0.004069 -2.050 0.040405 *
## inter.travel.red1 -0.030012 0.061666 -0.487 0.626478
## inter.travel.red2 -0.162898 0.067333 -2.419 0.015551 *
## wear.mask.red1 -0.238948 0.063963 -3.736 0.000187 ***
## wear.mask.red2 -0.136099 0.076365 -1.782 0.074714 .
## wear.mask.red3 -0.080739 0.060022 -1.345 0.178572
## test.policy.red1 -0.060338 0.057882 -1.042 0.297213
## test.policy.red2 -0.101022 0.068486 -1.475 0.140189
## time -0.034176 0.004858 -7.036 1.99e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
p_load(plotly, modelbased)
data0.anal1$Predicted <- estimate_response(model.glm)$Predicted
## Warning in get_predicted.glmmTMB(model, data = data, predict = predict, : `predict = 'prediction'` is currently not available for glmmTMB models.
## Changing to `predict = 'expectation'`.
plot2 <- data0.anal1 %>%
ggplot() +
geom_line(aes(x = log(week.case.eff), y = log(week.case.eff)), linetype = "dashed") +
geom_point(aes(x = log(week.case.eff) , y = log(Predicted), key=iso_code), color = "red") +
ylab("wADGR (predicted)") + xlab("wADGR (observed)")
## Warning: Ignoring unknown aesthetics: key
plot2
ggplotly(plot2, source = "select", tooltip = c("key") )
# Develop function marginal
amex<- function(ame){
x1 <- data.frame(list("AME" = ame*100, "factors" = as.character(substitute(ame))))
}
summary(model.glm)
## Family: beta ( probit )
## Formula:
## week.case.eff ~ uhc_index + inter.travel.red + wear.mask.red +
## test.policy.red + time + (time | iso_code)
## Data: data0.anal
##
## AIC BIC logLik deviance df.resid
## -2679.9 -2619.9 1354.0 -2707.9 536
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev. Corr
## iso_code (Intercept) 0.1224981 0.3500
## time 0.0002592 0.0161 -0.79
## Number of obs: 540, groups: iso_code, 30
##
## Dispersion parameter for beta family (): 43.3
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.654016 0.282376 -2.316 0.020552 *
## uhc_index -0.008340 0.004069 -2.050 0.040405 *
## inter.travel.red1 -0.030012 0.061666 -0.487 0.626478
## inter.travel.red2 -0.162898 0.067333 -2.419 0.015551 *
## wear.mask.red1 -0.238948 0.063963 -3.736 0.000187 ***
## wear.mask.red2 -0.136099 0.076365 -1.782 0.074714 .
## wear.mask.red3 -0.080739 0.060022 -1.345 0.178572
## test.policy.red1 -0.060338 0.057882 -1.042 0.297213
## test.policy.red2 -0.101022 0.068486 -1.475 0.140189
## time -0.034176 0.004858 -7.036 1.99e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
coef_model <- summary(model.glm)$coefficients$cond[,1][-1]
coef <- as.data.frame(coef_model)
# AME time
AME.time <- coef_model[9] * mean(dnorm(predict(model.glm, type="link")))
z9 <- amex(AME.time)
# uhc index
AME.uhc_idx <- coef_model[1]* mean(dnorm(predict(model.glm, type="link")))
z1 <- amex(AME.uhc_idx)
# travel
trv0 <- predict(model.glm, newdata=data.frame(list("time"=data0.anal1$time,
"uhc_index"=data0.anal1$uhc_index,
"inter.travel.red"="0",
"wear.mask.red"=data0.anal1$wear.mask.red,
"test.policy.red"=data0.anal1$test.policy.red,
"iso_code"=data0.anal1$iso_code) ), type = "response")
trv1 <- predict(model.glm, newdata=data.frame(list("time"=data0.anal1$time,
"uhc_index"=data0.anal1$uhc_index,
"inter.travel.red"="1",
"wear.mask.red"=data0.anal1$wear.mask.red,
"test.policy.red"=data0.anal1$test.policy.red,
"iso_code"=data0.anal1$iso_code) ), type = "response")
trv2 <- predict(model.glm, newdata=data.frame(list("time"=data0.anal1$time,
"uhc_index"=data0.anal1$uhc_index,
"inter.travel.red"="2",
"wear.mask.red"=data0.anal1$wear.mask.red,
"test.policy.red"=data0.anal1$test.policy.red,
"iso_code"=data0.anal1$iso_code) ), type = "response")
trv10 <- trv1 - trv0
trv20 <- trv2 - trv0
AME.trv10<-mean(trv10)
AME.trv20<-mean(trv20)
z2<-amex(AME.trv10)
z3<-amex(AME.trv20)
#wear.mask.red
mask0 <- predict(model.glm, newdata=data.frame(list("time"=data0.anal1$time,
"uhc_index"=data0.anal1$uhc_index,
"inter.travel.red"=data0.anal1$inter.travel.red,
"wear.mask.red"="0",
"test.policy.red"=data0.anal1$test.policy.red,
"iso_code"=data0.anal1$iso_code) ), type = "response")
mask1 <- predict(model.glm, newdata=data.frame(list("time"=data0.anal1$time,
"uhc_index"=data0.anal1$uhc_index,
"inter.travel.red"=data0.anal1$inter.travel.red,
"wear.mask.red"="1",
"test.policy.red"=data0.anal1$test.policy.red,
"iso_code"=data0.anal1$iso_code) ), type = "response")
mask2 <- predict(model.glm, newdata=data.frame(list("time"=data0.anal1$time,
"uhc_index"=data0.anal1$uhc_index,
"inter.travel.red"=data0.anal1$inter.travel.red,
"wear.mask.red"="2",
"test.policy.red"=data0.anal1$test.policy.red,
"iso_code"=data0.anal1$iso_code) ), type = "response")
mask3 <- predict(model.glm, newdata=data.frame(list("time"=data0.anal1$time,
"uhc_index"=data0.anal1$uhc_index,
"inter.travel.red"=data0.anal1$inter.travel.red,
"wear.mask.red"="3",
"test.policy.red"=data0.anal1$test.policy.red,
"iso_code"=data0.anal1$iso_code) ), type = "response")
mask10<-mask1-mask0
mask20<-mask2-mask0
mask30<-mask3-mask0
AME.mask10<-mean(mask10)
AME.mask20<-mean(mask20)
AME.mask30<-mean(mask30)
z4<-amex(AME.mask10)
z5<-amex(AME.mask20)
z6<-amex(AME.mask30)
# test
test0<- predict(model.glm, newdata=data.frame(list("time"=data0.anal1$time,
"uhc_index"=data0.anal1$uhc_index,
"inter.travel.red"=data0.anal1$inter.travel.red,
"wear.mask.red"=data0.anal1$wear.mask.red,
"test.policy.red"="0",
"iso_code"=data0.anal1$iso_code) ), type = "response")
test1<- predict(model.glm, newdata=data.frame(list("time"=data0.anal1$time,
"uhc_index"=data0.anal1$uhc_index,
"inter.travel.red"=data0.anal1$inter.travel.red,
"wear.mask.red"=data0.anal1$wear.mask.red,
"test.policy.red"="1",
"iso_code"=data0.anal1$iso_code) ), type = "response")
test2<- predict(model.glm, newdata=data.frame(list("time"=data0.anal1$time,
"uhc_index"=data0.anal1$uhc_index,
"inter.travel.red"=data0.anal1$inter.travel.red,
"wear.mask.red"=data0.anal1$wear.mask.red,
"test.policy.red"="2",
"iso_code"=data0.anal1$iso_code) ), type = "response")
test10<-test1-test0
test20<-test2-test0
AME.test10<-mean(test10)
AME.test20<-mean(test20)
z7<-amex(AME.test10)
z8<-amex(AME.test20)
z_merge <- rbind(z1,z2,z3,z4,z5,z6,z7,z8,z9)
term3 <- c(
"UHC Index\n",
"Ban arrivals from\nsome regions\n",
"Ban on all regions\n",
'Recommend or Required wear\nmask in some public spaces',
'Required wear mask in all\n public spaces with \nother people present',
'Required wear mask outside\n the home at all times',
"Testing of anyone \nshowing symptoms\n",
"Open public testing\n",
"Time\n"
)
p_load(broom.mixed)
sum_table <- tidy(model.mul, conf.int = T)
# create new table summary of model
new_table <- sum_table[2:10,]
new_table1 <- cbind(new_table, term3, z_merge)
new_table2 <- new_table1 %>% select(term, term3, AME, estimate, conf.low, conf.high)
new_table2$term3 <- factor(new_table2$term3, levels = new_table2$term3)
p.vac.period<- ggplot(new_table2, aes(x = term3, y = estimate,
ymin = conf.low, ymax = conf.high)) +
geom_hline( yintercept = 0, color = 'Black' ) + ylim(-0.5,0.5) +
geom_linerange(size =1, color = 'steelblue') + geom_point(size=1.5, color ='steelblue') + coord_flip() +
ylab("Probit Transformation of Average Weekly Growth Rate ") +
xlab("Predictors") +
theme_minimal() +
theme(strip.background = element_rect(NA),
axis.text = element_text(size = 20),
axis.title = element_text(size = 20),
legend.text = element_text(size = 10),
legend.title = element_text(size = 10),
strip.text = element_text(size = 8),
legend.position = "top")
p.vac.period
ame.plot <- ggplot(data = new_table2, aes(x= term3, y = AME)) + geom_bar(stat="identity",fill="steelblue") + coord_flip() +
geom_hline( yintercept = 0, color = 'Black' ) +
geom_text(aes(label=round(AME,2)), hjust= 1.05, color="black", size=5)+
ylab("Average marginal effect (%)") +
ylim(-3,3) +
xlab("") +
theme_minimal() +
theme(strip.background = element_rect(NA),
axis.text.y=element_blank(),
axis.text = element_text(size = 20),
axis.title = element_text(size = 20),
legend.text = element_text(size = 10),
legend.title = element_text(size = 10),
strip.text = element_text(size = 8),
legend.position = "top")
ame.plot
cairo_pdf(filename = "F:/STUDY HCA IN TAIWAN/0. ideas for thesis/COVID19 and NPIs/Data analysis/results/plot_sep/fig3_2020_Sep29.pdf",
width = 20, height = 10, #inch
onefile = TRUE, family = "Segoe UI")
merge_bar <- ggarrange(p.vac.period, ame.plot, ncol = 2, nrow = 1, widths = c(1.5,1))
merge_bar
dev.off()
## png
## 2